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PREFACE 


NASA/Marshall Space Flight Center (MSFC) has developed a new atmospheric turbulence 
model which is more realistic and less conservative when applied in space shuttle reentry simula- 
tions involving engineering calculations of reaction control system fuel expenditures. Both Georgia 
Institute of Technology and BDM Corporation were contracted to update the required turbulence 
velocity sigmas and length scales, and to apply them in a white noise filter technique to arrive at a 
more realistic engineering turbulence model. This model is also envisioned to be useful in other 
type spacecraft and aircraft simulation studies. This project was funded by the NASA-Johnson 
Space Center Space Shuttle Office. 
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TECHNICAL MEMORANDUM 


NEW ATMOSPHERIC TURBULENCE MODEL 
FOR SHUTTLE APPLICATIONS 

I. INTRODUCTION 


It has been determined at NASA [1] that the currently used atmospheric turbulence (wind) 
model [2,3] for space shuttle reentry simulation is overly conservative. Use of this model in shuttle 
reaction engine fuel usage calculations assumes severe turbulence all the way from reentry to land- 
ing. Johnson Space Center (JSC) sets reaction control system redlines based on these fuel predic- 
tions. However, in reality, the orbiter generally returns from space with approximately 270 kg (600 
lb) of extra, unused fuel aboard creating an unneeded weight excess. Turbulence in the real 
atmosphere is patchy or intermittent with quiescent zones present. Therefore, the Environmental 
Analysis Branch of Marshall Space Flight Center’s (MSFC) Earth Science and Applications Divi- 
sion developed a more realistic engineering-oriented turbulence model that can now be used in 
shuttle simulation work to select a more rational fuel reserve redline value, along with other poten- 
tial atmospheric turbulence applications. 

This modeling task was accomplished in two parts. First, Dr. C.G. Justus of the Georgia 
Institute of Technology updated the statistical turbulence data base by a literature search to arrive at 
better estimates of anisotropic horizontal and vertical turbulence velocity standard deviations (a u 
and cr w ), and length scale parameters (L x and L z ), from near-surface to 200-km altitude. The 
y-component (v) was not explicitly calculated, but is generated identically to the z-component. 

These model statistics are available in the form of a program subroutine to evaluate turbulence a’s 
and L’s as a function of altitude. This task is fully documented as part II of this report, and was 
taken completely from the Georgia Tech final report [4] and reprinted in its entirety here. The 
results from this task have also been presented at a recent technical conference [5]. 

The second part of this task was done by Dr. C.W. Campbell and M.K. Doubleday of 
BDM Corporation who applied and modeled the new turbulence statistics of Justus [4] in a proce- 
dure which inputs Gaussian white noise into a low-pass linear filter to output the simulated turbul- 
ence in a Gaussian time series. The transfer function of the filter was selected to produce a desired 
von Karman spectrum, with a more realistic probability distribution. In Campbell’s study, for long- 
itudinal spectra, transfer function approximations to the von Karman transfer function, up to fifth 
order, were derived versus differing sampling rates. The corresponding transverse transfer function 
is one order higher. The resulting longitudinal and transverse equations used can be directly coded 
into the shuttle reentry simulation, or into any other type of vehicle flight simulation procedure. 

The only inputs to the equations are the appropriate turbulence length scale and relative wind 
velocity turbulent intensity and the sampling rate. Campbell’s work is documented in the BDM 
final report [6] which is also presented in this report as part III. This work has also been presented 
at technical conferences [7,8]. Larry McWhorter at NASA/JSC is currently implementing this new 
turbulence model in his shuttle reentry fuel simulation work. 
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II. UPDATED TURBULENCE STATISTICS AND SUBROUTINE 


A. BACKGROUND 


As evidenced by the information in Table 1 , from a review of turbulence models by 
Moorhouse and Heffley (1986), there are a wide variety of techniques which have been used in 
turbulence simulation. The purpose of this study was to develop a turbulence model for use in 
space shuttle reentry simulations, which will be simple to use, computationally fast, and consistent 
with techniques of Monte Carlo modeling and digital filtering, being developed by other researchers 
for NASA (Campbell and Fichtl, 1985; Campbell, 1986). 

Previous NASA methods for simulation of turbulence for aerospace vehicle flight simula- 
tions include the turbulence simulation technique of Fichtl (1977) and the turbulence criteria and 
model presented in Section 2.4 of Turner and Hill (1982). The Fichtl approach was used in 
Appendix 10.10 of “Natural Environment Design Requirements for the Space Shuttle” (NASA, 
1975), and formed the basis for the shuttle simulation turbulence tapes of Tatom and Smith (1982). 
The Turner and Hill approach has also been adopted and recommended as the turbulence design 
criteria for other NASA projects (e.g., Adelfang, 1987). 

The turbulence model proposed here is based on these techniques, with updates and 
modifications. The method of Turner and Hill is based on a model probability distribution p(a) 
given by: 


p(o-) = V2 /tt (P,/bj) exp 



+ V2/tt (P 2 /b 2 ) exp 



( 1 ) 


(their equation 2.38), where b) and b 2 are the standard deviations of rms gust velocity in nonstorm 
and storm turbulence, and Pi and P 2 are the fractions of flight time or distance flown in nonstorm 
and storm turbulence. Equation (1) assumes a fraction P 0 for flight time or distance in smooth air, 
such that 


P 0 + P, + P 2 = 1 


(2) 


It should be noted that, for consistency with equation (2), equation (1) should have an added term 
P 0 5(or), where 8 is the Dirac delta function. 
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B. BASIS OF THE REVISED MODEL 


There are several changes in the form of the model developed. The half-Gaussian distribu- 
tion of equation (1) has the unrealistic feature that the most probable value (the mode) of cr in the 
p(a) distribution is a = 0. The alternate form of distribution suggested is the Rayleigh, which has 
a more probable value (mode) closer to the mean (expected) value of the distribution. The Rayleigh 
distribution is given by 


p(or) = (cr/b 2 ) exp [-(a/b) 2 /2] 


(3) 


which has a mean value cr = (ir/2) l/2 b, a mode of b, and a standard deviation of [2 - ir/2] l/2 b. 

For the cumulative probability p(a^a,), the Rayleigh distribution produces a Gaussian distribution 

p(o-^ai) = exp[-(o-|/b) 2 /2] . (4) 


For implementation in the Monte Carlo series simulation (Campbell and Fichtl, 1985; 
Campbell, 1986), the values of cr can be evaluated from Gaussian-distributed components, <j\ and 
cr j (not to be confused with spatial vector components) by the relation cr = (cr { 2 + <jj 2 ) l/2 , where 
CTj and a, are each selected from a Gaussian (normal) distribution, that is 


p(o-j) = (2irb 2 ) 1/2 exp[-(oyb) 2 /2] 


(5) 


and similarly for p(cxj). 

Equation (3) is used separately to evaluate the probability distribution for no turbulence 
(with cr = 0 and probability P 0 ), moderate turbulence (with cr selected from the distribution with b 
= b| and probability P|), and severe turbulence (with cr selected from the distribution with b = b 2 
and probability P 2 ). Instead of using these probability distributions additively, as in equation (1), 
each of these distributions is used separately, according to which severity level of turbulence is 
being encountered (none, moderate, or severe). A turbulence severity parameter a (with value 
selected from a uniform random distribution between 0 and 1) is used to determine the severity 
level of the turbulence: there is no turbulence if a < P 0 ; the turbulence is severe if a s? 1 - P 2 , 
and moderate otherwise. Minimum vertical depths for layers of moderate and severe turbulence are 
also specified. Thus, once the series simulation enters a zone with moderate or severe turbulence, 
it must remain at this severity level until at least the specified minimum depth (or the specified 
minimum horizontal extent) has been traversed, before it can return to a lower severity level. 


4 



The effects of vertical correlations between a values at one altitude with those at adjacent 
altitudes are also included. This feature incorporates the fact that rms turbulence gust magnitudes cr 
must vary more-or-less continuously along the trajectory [not discontinuously as if selected by 
independent calls upon the probability distribution of equation (3)]. The turbulence cr value changes 
abruptly, however, when transitioning from one severity layer to another (with only a one-step, 
linear interpolation smoothing being applied each time a new intensity layer is encountered). No 
correlations are assumed between sigmas in two layers which are spatially separated by a layer of 
lower turbulence magnitude. Thus, each time a layer of higher than current turbulence severity is 
encountered, the random number generator sequence for the sigma selections is reinitialized. 

A major portion of the model development project has been a literature survey to develop 
revised parameter values for the data on turbulence intensities (mean a values), scales, and proba- 
bilities of intensity levels. In addition, new parameter values were required for the vertical scale of 
the a interlevel correlation, and the minimum vertical sizes for moderate and severe turbulence 
layers. Anisotropic horizontal and vertical values are provided for the turbulence intensity and scale 
parameters. 


C. MODEL DESCRIPTION 


The literature survey for updating the turbulence parameter values consisted of a search of 
the Scientific and Technical Aerospace Abstracts, the International Aerospace Abstracts, and the 
Meteorological and Geoastrophysical Abstracts for the period 1970 to present. Parameter values 
resulting from this literature review have been incorporated into the turbulence simulation 
subroutine. 

Information on horizontal and vertical scales of turbulence (integral scale in the velocity 
correlation function) were averaged from data in McCloskey et al. (1971), U.S. Department of 
Defense (1975), Fichtl (1977), Hasty (1977), Justus et al. (1980), Turner and Hill (1982), Frost et 
al. (1985), Murrow (1986), and Reid and Vincent (1987). Resulting average values for these hori- 
zontal and vertical scales are shown in Figure 1 as a function of altitude from the surface to 200 
km. 


Values of the sigma components crj and oj to be selected from the Rayleigh distribution by 
equation (4) are assumed to be correlated over horizontal and vertical separations with scale values 
(integral scale) which are related to the horizontal and vertical scales of the turbulent velocity, 
through a ratio which was assumed to vary from 10.0 at the surface to 5.0 at 20 km and higher. 
For example, with a horizontal turbulence velocity scale of 0.52 km at the surface, the horizontal 
sigma scale is assumed to be 5.2 km. Assumptions such as these were required because no obser- 
vational data were found to provide direct estimates of the sigma scales. These ratio values (5-10) 
were assumed since turbulent velocity statistics cannot accurately be determined from measurements 
unless the turbulence is stationary (relatively constant sigma) over at least 5-10 velocity scale 
values; hence a presumption that the sigma scales relative to the velocity scales are of at least this 
magnitude range (with larger ratio values anticipated to occur at lower altitudes). 
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For data on the average value d u (the average longitudinal turbulence intensity), data were 
averaged over information taken from Pershikov (1969), McCloskey et al. (1971), Ryan et al. 
(1971), Fichtl (1977), Kao et al. (1977), Waco (1978), Justus et al. (1980), Vinnichenko et al. 
(1980), Moorhouse and Woodcock (1982), Murrow et al. (1982), Turner and Hill (1982), Hocking 
(1983), Frost et al. (1985), Hill (1986), Murrow (1986), and Andrews et al. (1987). For the Fichtl 
(1977) values of ct u , a weight value of 0.6 was used, since these values were taken to represent 
extreme values rather than average values. 

Data on the average value <r w (or on the ratio of cr w /a u ) were taken from McCloskey et al. 
(1971), Ryan et al. (1971), U.S. Department of Defense (1975), Fichtl (1977), Waco (1978), 
Vinnichenko et al. (1980), Murrow et al. (1982), Turner and Hill (1982), Frost et al. (1985), Sato 
(1985), and Murrow (1986). All of these data were interpreted in terms of an average ratio for a w / 
a u , which was then multiplied by the average tx u , determined above, to arrive at the final value for 
the average ct w . 

Values of the ratio (y = ct 2 /ct 0 of the magnitude of severe turbulence (a 2 ), relative to the 
magnitude of moderate turbulence (<T|), were averaged from data in Kao et al. (1977), Moorhouse 
and Woodcock (1982), and Turner and Hill (1982). No values for the y intensity ratio were found 
for altitudes above 18 km; therefore a value of <r 2 equal to 1.5 times <j t was assumed at 30 km 
and higher, with a smooth transition from the value of y = 1.84 at 18 km to y = 1.5 at 30 km. 

Values of the total probability P = Pj + P 2 for encountering turbulence (either moderate or 
severe) were taken and averaged from McCloskey et al. (1971), Wilson et al. (1971), Ehrenberger 
(1975), Waco (1976), Hasty (1977), Zimmerman and Murphy (1977), Turner and Hill (1982), and 
Ehrenberger (1987). Values for p = P 2 /P|, the ratio of the probability of encountering severe 
turbulence to that for encountering moderate turbulence, were obtained and averaged from 
McCloskey et al. (1971), Waco (1976), and Turner and Hill (1982). No values of p were found in 
the literature for altitudes above 30 km (where p = 0.1). A steady decrease to a value of p = 

0.01 was assumed at 120 km, with p = 0.01 between 120 and 200 km. From values of the total 
probability P, and the ratio p, separate values of P| and P 2 were obtained from 


Pi = P/(l + p) and P 2 = pP/(l + p) 


( 6 ) 


From values of p and the severity magnitude ratio y, separate values for average sigma for mod- 
erate turbulence ((X|) and for severe turbulence (ct 2 ) can be obtained from the overall average sigma 
a for both horizontal and vertical components (u and w) via 


a lu 


(1 + P)^u 
(1 + py) 


and 


°2u 


7(1 + p)q~u 
(1 + py) 


(7) 


^ 1 W 


(1 + pKv 
(1 + py) 


and 


y(l + p)<?w 
(1 + py) 


( 8 ) 
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Values for the resulting averages ct 2u and cr 2w as a function of altitude for severe turbulence 
are shown in Figure 2. The profiles of a lu and ct 1w values for moderate turbulence are shown in 
Figure 3. Values of the average sigmas for composite turbulence, cf c , defined to be probability- 
weighted values of cf| and cr 2 , are defined by 


Cc — Po'(0) + Pi or, + P 2 a 2 


(9) 


for both u and w components [where the normalization P 0 + Pi + P 2 = 1 has been used, from 
equation (2)]. The resulting composite turbulence cr c values are shown in Figure 4, as functions of 
height, for both horizontal and vertical components. 

Up to about 20 km, the magnitudes (Figures 2-4) and scales (Figure 1) for the turbulence 
model are taken to represent atmospheric turbulence as conventionally defined. Above this height, 
the magnitudes and scales are taken to represent either gravity waves (which in high-speed flight 
may produce the same dynamical effects on the vehicle as does turbulence), or the turbulence 
which may result from gravity wave breaking or gravity wave dissipation at critical levels. 

For the minimum vertical extent of turbulence layers, data values from Ehrenberger (1975), 
Rottger (1980), Barat (1982), Tanaka and Yamanaka (1984), Yamanaka and Tanaka (1984), Cot 
and Barat (1986), and Ehrenberger (1987) were averaged. No data on vertical layer thickness were 
available above the 60-90 km layer, so an extrapolated increase rate in the vertical layer size 
above this altitude range was assumed. 


Horizontal minimum extent of turbulence layers was derived from average data values taken 
from Ehrenberger (1975), Murrow et al. (1982), Hasty (1977), and Waco (1978). No data were 
available above 25 km, so an extrapolated height increase for the horizontal layer size was assumed 
above this height. No clear data were available on the relative sizes for severe turbulence relative 
to the size for moderate turbulence layers. Therefore the minimum vertical and horizontal sizes for 
severe turbulence layers were taken to be 1/2 that for the moderate layers. 

The attached program listing shows the model implementation as subroutine TURBSIG, 
which is to be used for simulating turbulence intensity values (cr’s) at a given altitude z (0 =s z «£ 
200 km). The subroutine returns values for cr u and cr w , the horizontal and vertical sigma 
components (in m/s), and L x and L z , the horizontal and vertical integral scale values (in km). 

Several model options (specified by the input parameter MODEL) are available: MODEL = 
0 for the complete model with all specified parameters (probabilities P 0 , Pi, P 2 , average sigma 
values, a lu , cr 2u , a lw , a 2w , scales L x , L z , etc.) with the local, spatially dependent, sigma values 
selected from the appropriate Rayleigh distribution for moderate and severe turbulence intensity 
levels; MODEL = 1 for moderate turbulence throughout the whole run, with sigma values selected 
from the Rayleigh distribution which has a = ff t ; MODEL = 2 for severe turbulence throughout 
the whole run with sigma values selected from the Rayleigh distribution which has a = a 2 ; 
MODEL = 3 for composite turbulence (probability-weighted average of severe, moderate, and 
non-turbulent) throughout the whole run, with sigmas selected from the Rayleigh distribution which 
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has ct = ct c ; MODEL = 4 for no turbulence throughout the whole run (cr = 0); MODEL = 5 for 
moderate turbulence with a = <X| (no Rayleigh distribution) throughout the whole run; MODEL = 
6 for severe turbulence with a = <x 2 throughout the whole run (no Rayleigh distribution); and 
MODEL = 7 for composite turbulence with cr = ct c throughout the whole run (no Rayleigh 
distribution). 

The main program TESTSIG is designed to be used for producing test output of the 
TURBSIG subroutine, and illustrates how TURBSIG is to be used in the reentry trajectory simula- 
tion programs. The user must select a starting random number (any odd positive integer) for the 
random number generator (RAND), and a value of MODEL must be selected. The TURBSIG 
subroutine must also be called once with displacement values dx = 0 and dz = 0, in order to 
initialize all of its variables. On subsequent calls, the trajectory program must pass to the subrou- 
tine the displacement values (dx and dz, in km) of current position from previous position. Only 
the magnitudes of the displacements are of importance, so positive or negative values may be used. 
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Severe Turbulence Magnitude 



Figure 2. Horizontal and vertical magnitudes for severe turbulence as a function of height 
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Figure 3. Horizontal and vertical magnitudes for moderate turbulence as a function of height. 
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Figure 4. Horizontal and vertical magnitudes for composite turbulence (probability-weighted severe, 
moderate, and no turbulence cases) as a function of height. 







TABLE 1 . A SURVEY OF ATMOSPHERIC DISTURBANCE MODELS 
(Moorhouse and Heffley, 1986) 


Model 

Key Features 

Sources 

Dryden turbulence 

A convenient spectral form based on an 
exponential autocorrelation function 
for the axial component 

Dryden (1961) 

von Karman turb- 
ulence 

A spectral form for which the correla- 
tion function includes a finite micro- 
scale, thus the relative proportion of 
spectral power at high frequencies ex- 
ceeds that of the Dryden. 

von Karman 
(1961), Hou- 
bolt (1973) 

Ornstein-Uhlenbeck 

turbulence 

A spectral form with first-order long- 
itudinal and transverse components 

Gaonkar 

(1980) 

Etkin one dimen- 
sional turbulence 
power spectra 

The local turbulent velocity field is 
approximated by a truncated Taylor 
series which yields uniform and gradi- 
ent components. High frequency spec- 
tral components eliminated on the 
basis of aircraft size. Based on Dry- 
den form, but gradient spectra are 
non-realizable unless simplified. 

Etkin (1961), 
Etkin (1959a, 
b), Etkin 
(1972) 

Versine gust 

A discrete gust waveform 

Moorhouse and 

Woodcock, 

(1982) 

Lappe 1 ow-al t i t ude 
turbulence model 

Experimentally-obtained data of vert- 
ical gust spectra, mean wind speed, 
and lapse rate were used to develop a 
low-level turbulence model. The turb- 
ulence spectra are presented for diff- 
erent types of terrain, height, and 
meteorological conditions. 

Lappe (1966) 

Multiple point 
source turbulence 

A two-dimensional gust field generated 
from two or more noise sources having 
prescribed correlation functions and 
located spanwise or lengthwise on the 
vehicle. 

Etkin (1980), 
Holley and 
Bryson (1975), 
Skelton (1968) 



TABLE 1. (Continued) 


Model 

Key Features 

Sources 

Holley-Bryson 
random turbulence 
shaping filters 

A matrix differential equation formu- 
lation of uniform and gradient compo- 
nents including aircraft size effects. 
Filter equation coefficients deter- 
mined from least squares fit to multi- 
point-source-derived correlation func- 
tions . 

Holley and 
Bryson (1975) 

University of 
Washington non- 
Gauss ian atmos- 
pheric turbulence 
model 

Non-Gaussian model using modified Bes- 
sel functions to simulate the patchy 
characteristics of real-world turbu- 
lence. Spectral properties are Dryden 
and include gust gradients. 

Reeves et al. 
(1974), Reeves 
(1969) 

Delft University 
of Technology non- 
Gauss ian structure 
of the simulated 
turbulent environ- 
ment 

Non-Gaussian model similar in form to 
the University of Washington model, 
but uses the Hilbert transformation to 
model intermittency as well as patchi- 
ness. Includes University of Washing- 
ton model features extended to approx- 
imate transverse turbulence velocities 
and gradients. 

van de Moes- 
dijk (1978) 

Royal Aeronautical 
Establishment 
model of non-Gaus- 
sian turbulence 

Non-Gaussian turbulence model with a 
variable probability distribution 
function and a novel digital filtering 
technique to simulate intermittency. 
Spectral form approximately von 
Karman. 

Tomlinson 
(1975), Jones 
(1976), Jewell 
and Heffley 
(1978) 

The Netherlands 
National Aerospace 
Laboratory model 
of non-Gaussian 
turbulence 

Similar to the Royal Aeronautical 
Establishment model, but extended to 
include patchiness and gust gradient 
components and transverse velocities. 

Jansen (1977a, 
1977b) 

University of 
Virginia turbu- 
lence model 

Models patchiness by randomizing gust 
variance and integral scale of basic 
Dryden turbulence. 

Jacobson and 
Joshi (1977) 
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TABLE 1. (Continued) 


Model 

Key Features 

Sources 

MIL Standard 
turbulence model 

First order difference equation imple- 
mentation of turbulence filters based 
on 8785 Dryden turbulence and refitted 
rolling gust intensity. 

Hoh et al. 
(1982) 

Indian Institute 
of Science non- 
stationary turbu- 
lence model 

Non-stationary turbulence is obtained 
over finite time windows by modulating 
a Gaussian process with either a de- 
terministic or random process. The 
result is patchy-like turbulence, sim- 
ilar to the University of Washington 
model, except the time- varying statis- 
tics of turbulence are presented for 
the deterministic modulating func- 
tions. 

Gaonkar (1980) 

FAA wind shear 
models 

Three-dimensional wind profiles for 
several weather system types including 
fronts, thunderstorms, and boundary 
layer. The profiles are available in 
table form. 

Foy and Gart- 
ner (1979), 
Frost and Camp 
(1977) 

STI wind shear 
model 

Time and space domain models of mean 
wind and wind shear (ramp wave forms) 
are combined with MIL-F-8785C Dryden 
turbulence to obtain the total atmos- 
pheric disturbance. The magnitudes of 
the mean wind and wind shear are eval- 
uated in terms of the aircraft’s ac- 
celeration capabilities. 

Hoh and Jewell 
(1976), Hef- 
fley and Jew- 
ell (1978) 

Sinclair frontal 
wind shear model 

A generic model of frontal surface 
wind shear derived from a reduced- 
order of Navier-Stokes equations. 
Relatively simple to use and can match 
the overall characteristics of meas- 
ured wind shears. 

Jewell et al. 
(1979), Sin- 
clair and West 
(1978) 


21 




TABLE 1. (Continued) 


Model 

Key Features 

Sources 

MIL-F-8785B atmos- 
pheric disturbance 
model 

Intensities and scale lengths are 
functions of altitude and use either 
Dryden or von Kannan spectral forms 
or a one minus cosine gust. Also 
spectral descriptions of rotary gusts. 

U.S. Dept, of 
Defense 
(1969), Chalk 
et al. (1969) 

MIL-F-8785C atmos- 
pheric disturbance 
model 

Same as 8785B with the addition of a 
logarithmic planetary boundary layer, 
a vector shear and a Naval carrier 
airwake model. 

U.S. Dept, of 
Defense (1980) 

ESDU atmospheric 
turbulence 

Rather general, but contains compre- 
hensive descriptive data for turbu- 
lence intensity, spectra, and proba- 
bility density. 

Anonymous 
(1974, 1975) 

Boeing atmospheric 
disturbance model 
turbulence 

A comprehensive model of atmospheric 
disturbances that includes mean wind, 
wind shear, and random turbulence. 
Turbulence is Gaussian and uses linear 
filters that closely approximate the 
von Karman spectral form. Mean wind 
and turbulence intensity are functions 
of meteorological parameters. 

Barr et al. 
(1974) 

Wasicko carrier 
airwake model 

Includes mean wind profile, effect of 
ship motion, and turbulence. 

Durand (1967) 

Naval ship airwake 
model 

Includes free air turbulence filters 
plus steady, periodic, and random com- 
ponents of airwake which are functions 
of space and time. 

U.S. Dept. 
Defense 
(1980), Nave 
(1978) 

Vought airwake 
model for DD-963 
class ships 

Combined random and deterministic wind 
components for free air and ship air- 
wake regions. Based on wind tunnel 
flow measurements. 

Fortenbaugh 

(1978) 
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TABLE 1. (Concluded) 


Model 

Key Features 

Sources 

STI wake vortex 
encounter model 

A two-dimensional model of the flow 
field due to the wake vortex of an 
aircraft is presented. The parameters 
of the flow-field model are weight, 
size, and speed of the vortex-generat- 
ing aircraft, and distance and orien- 
tation of the vortex-encountering 
aircraft. Strip theory is used to 
model the aerodynamics of the vortex- 
encountering aircraft. 

Johnson and 
Teper (1974) 

Campbell and 
Sandborne wind 
shear and turbu- 
lence model 

Spatial model based on joint airport 
weather studies (JAWS) microburst 
data. Permits calculation of aero- 
dynamic loads over body of aircraft. j 

Campbell and 

Sandborne 

(1985) 

Zhu and Etkin 
microburst model 

Generic spatial model of microburst 
velocity components based on poten- 
tial flow singularity distribution 
involving only three adjustable para- 
meters. 

Zhu and Etkin 
(1985) 
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LISTING FOR TESTSIG PROGRAM TO TEST TURBULENCE SIMULATION MODEL FOR 
SIGMA VALUES AND TURBSIG SIMULATION MODEL SUBROUTINE 


C PROGRAM - TESTSIG to test the TURBSIG subroutine 

implicit double precision (a-h,o~z) 

Open (21, File- ’sigmau’ , Status=’New’ ) 

Open (22, File= ’ sigmaw’ , Status=’New* ) 

Open(23,File=’scaleu’ , Status= ’ New’ ) 

Open(24, File= ’ scalew* , Status- * New* ) 

Open(25,File=’ severity’ , Status^ ’New’ ) 

Open (26, Fi le= ’model ’ , Status^’New’ ) 

1 Write(*,5) 

5 Format (’ Enter starting random number (odd integer): ’ ) 
Read(*, *)nl 
if (nl. le. 0)goto 99 
if (MOD(nl, 2) . ne. l)goto 1 
z = rand(nl) 

Write(*, 10) 

10 Format (’ Enter model number 0-7: ’ ) 

Read(*,*)model 
Write(*, 15) 

15 Format (’ Enter starting height: * ) 

Read(*,*)Z 

Call turbsig(model,z,0. ,0. , sigmau, sigmaw, xlwind, zlwind, isev) 
Write(*, 30 )z, sigmau, sigmaw, xlwind, zlwind, isev 
Write(*, 20) 

20 Format ( 5 Enter displacements DX, DZ in km: ’) 

Read ( * , * ) dx , dz 
x = 0. 

xmax = 99999. 
if (abs(dz) . le. 0. ) then 
Write(*, 23) 

23 Format (’ Enter maximum x value to simulate: ’ ) 

Read (*, *) xmax 
endif 

25 Call turbsig(model , z , dx , dz , sigmau, sigmaw, xlwind, zlwind, isev) 
Write(*, 30) z , sigmau, sigmaw, xlwind, zlwind, isev 
Write (21 ,40) z, sigmau 
Write (22, 40) z , sigmaw 
Write ( 23 , 40 ) z , xlwind 
Write(24,40)z, zlwind 
Write(25,50)z, isev 

Wr ite( 26, 30) z, sigmau, sigmaw, xlwind, zlwind, isev 
30 Format (5f 10. 2, i3) 

40 Format (2f 10. 2) 

50 Format (f 10. 2, i4) 

Z = Z - abs(dz) 
x = x + abs(dx) 

if (Z.lt.O.O.or.x.gt.xmax)goto 99 
goto 25 
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99 stop 
end 


SUBROUTINE TURBSIG(MODEL,Z,DX,DZ,SIGMAU,SIGMAW,XLWIND,ZLWIND, 
& ISEV) 


c 

c. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Simulation of turbulent wind standard deviation, sigma, selected 
from a Rayleigh distribution, with parameters which are a 
function of height Z. 

Three turbulence intensities are simulated: a non~turbulent 
background, moderate turbulence, and severe turbulence. 
Frequencies of occurrence and minimum persistence of 
layers are specified for moderate and severe intensities. 
Correlation lengths for the turbulent wind field and for 
the sigmas for the wind field are also specified 
independently. 

The characteristics of the Rayleigh distribution 

for standard deviation (sigma) are determined by the expected 

value (average value) of the sigma distribution. 


Input subroutine arguments are: 

MODEL - 0 for complete model with all specified parameters; 
sigmas selected from Rayleigh distribution 

1 for moderate turbulence throughout the whole run; 

sigmas selected from Rayleigh distribution 

2 for severe turbulence throughout the whole run; 

sigmas selected from Rayleigh distribution 

3 for composite turbulence (average of severe, 
moderate and non turbulence) throughout the 
whole run; sigmas selected from Rayleigh 
distribution 

4 for no turbulence throughout the whole run 

5 for moderate turbulence of average sigma 
throughout the whole run 

6 for severe turbulence of average sigma 
throughout the whole run 

7 for composite turbulence of average sigma 
throughout the whole run 

Z ~ Current altitude in km 

DX - Horizontal displacement since last position, in km 
DZ - Vertical displacement since last position, in km. 

Note - To initialize values call TURBSIG with both 
DX = 0 and DZ = 0 
Output subroutine arguments are: 

SIGMAU - Current turbulence standard deviation for the 
horizontal wind components, in m/s 
SIGMAW - Current turbulence standard deviation for the 
vertical wind component, in m/s 
XLWIND - Current horizontal scale for turbulent wind, in km 
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ZLWIND - Current vertical scale for turbulent wind, in km 
ISEV - Severity parameter (1 = non turbulent, 

2 = moderate turbulence, 3 = severe turbulence, 
or 0 = composite turbulence: weighted average 
of categories 1-3) 


implicit double precision (a-h,o-z) 

The following static variables have values which should 
remain unchanged from one call of the subroutine to another. 
Some compilers require declaration of such static variables; 
others do not. 

Static sigmxu, sigmyu, sigmxw, sigmyw, sigsxu, sigsyu, 

& sigsxw, sigsyw, jsev, zsev, xsev, psigu, psigw, nsev 

double precision LSIGX(32) , LSIGZ(32) , LWINDX(32) , LWINDZ(32) , 

& LMODZ(32) , LSEVZ(32) , LMODX(32) , LSEVX(32) , maxz(3) ,maxx(3) 
dimension hgt ( 32 ) , S IGMXBAR ( 32 ) , S IGSXB AR( 32 ) , 

& SIGMZBAR(32) ,SIGSZBAR(32) ,PM(32) ,PS(32) ,sigu(3) ,sigw(3) 
data pi, one, two/3. 14159265359d0, 1 . OdO, 2. OdO/ 
data AFAC , BFAC/19. 51615854016301d0 , 1 . 00041693941245578d0/ 
Heights for turbulence model parameters, in km 
data hgt/0. , 1 . , 2. ,4. , 6. , 8. , 10. , 12 . , 14. , 16. , 18. ,20. ,25. , 30. , 
& 35. ,40. ,45. ,50. ,55. ,60. ,65. ,70. ,75. ,80. ,85. ,90. ,100. , 

& 120., 140., 160., 180., 200./ 

Mean Value for moderate turbulence sigmas, m/s (horizontal) 
Data SIGMXBAR/1. 25, 1.65,1.65,2.04,2.13,2.15,2.23,2.47, 

& 2.62,2.44,2.21,2.26,2.71,3.73,4.59,5.26,6.22,7.27,8.7, 

& 10.1,11.3,15.9,19.2,22.6,27.3,33.2,35.6,42.3,44.3,48.2, 

& 48.9,49.5/ 

Mean Value for moderate turbulence sigmas, m/s (vertical) 
Data SIGMZBAR/. 98, 1.36, 1.43, 1.68, 1.69, 1.69, 1.73, 1.79, 

& 1.91,2.10,2.07,1.99,2.09,2.39,2.58,2.87,3.25,4.21,4.40, 

& 4.42,4.05,5.04,6.3,8.3,10.3,11.8,11.4,10.7,10.8,11.7, 

& 11.8,12.0/ 

Mean Value for severe turbulence sigmas, m/s (horizontal) 
Data S IGSXB AR/3 .06,3.90,4.35,6.24,7.16,7.59,7.72,7.89, 

& 6.93,5.00,4.07,3.85,4.34,5.60,6.89,7.89,9.33,10.90,13.06, 

& 15.1,16.9,23.8,28.7,33.8,40.9,49.8,53.3,63.4,66.4,72.2, 

& 73.3,74.2/ 

Mean Value for severe turbulence sigmas, m/s (vertical) 

Data S IGSZB AR/2 .41,3.21,3.78,5.13,5.69,5.98,6.00,5.71, 

& 5.05,4.31,3.81,3.38,3.34,3.59,3.87,4.30,4.88,6.31,6.60, 

& 6.63,6.0,7.5,9.5,12.4,15.4,17.7,17.1,16.0,16.1,17.6, 

& 17.8,18.1/ 

Horizontal scale for turbulence sigmas, in km 

Data LSIGX/5. 2, 8. 3, 8. 6, 9. 4, 8. 8, 8. 3, 9. 2, 12. 6, 18., 21., 28., 

& 43. ,60. , 143. , 177. ,213. ,250. ,290. ,330. ,372. ,416. ,462. , 

& 5 10 . , 555 . , 605 . , 660 . , 765 . , 1 000 . , 1 160 . , 1350 . , 1500 . , 1500 . / 

Vertical scale for turbulence sigmas, in km 
Data LSIGZ/3. 2, 6. 2, 7. 7, 8. 5, 8. 4, 8. 1,8. 3, 11., 14., 16., 18. , 

& 22. ,33. ,44. ,42. ,31. ,26. ,27. ,30. ,34. ,38. ,41. ,45. ,49. , 

& 52. ,56. ,64. ,79. ,88. ,100. ,111. ,122./ 

Horizontal scale for turbulence winds, in km 

Data LWINDX/0. 520, 0.832, 0.902, 1.04, 1.04, 1.04, 1.23, 1.80, 
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& 2.82,3.40,5.00,8,64, 12.0,28,6,35.4,42.6,50. 1,57.9,66.0, 

& 74.4,83.2,92.3,102. ,111. ,121. ,132., 153. ,200. ,232. ,270. , 

& 300., 300./ 

Vertical scale for turbulence winds, in km 
Data LWINDZ/0. 323, 0.624, 0.831, 0.972, 1.01, 0.98, 1.10, 1.54, 

& 2.12,2.60,3.34,4.41,6.56,8.88,8.33,6.2,5.2,5.3,6.0,6.8, 

& 7.5,8.2,9.0,9.7,10.4,11.2,12.7,15.8,17.6,20.0,22.2,24.3/ 

c... Probability for encountering moderate turbulence 

Data PM/ . 867 ,. 199 ,. 0979 ,. 0738 , . 0650 , . 0704 , . 0677 , . 0502 , . 0368 , 
& . 0337 , . 0277 , . 0180 , . 0146, . 0185 , . 0249, . 0318, . 0386, . 0455 , 

& . 0682 ,. 0917 ,. 1620 ,. 2336 , . 3066 ,. 3810 ,. 5769 , . 7767 , . 9804 , 

& .9901,. 9901,. 9901,. 9901,. 9901/ 

c... Probability for encountering severe turbulence 

Data PS/ . 0 10 , . 025 , . 01 1 1 , . 0063 , . 0056 , . 0049 , . 0043 , . 0034 , 

& .0027, .0024, .0020, .0016, .0015, .0018, .0025, .0032, .0039, 

& . 0045 , . 0068, . 0083 ,. 0130 ,. 0164, . 0184, . 0190 , . 0231 , . 0233 , 

& .0196,. 0099 , . 0099 , . 0099 , . 0099 , . 0099/ 

c... Minimum vertical size for moderate turbulence layer, km 
Data LMODZ/0. 4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.64, 0.64, 0.64, 

& 0.64,0.50,0.40,0.56,0.72,0.88,1.04,1.20,1.36,1.48,1.64, 

& 1.80,1.96,2.12,2.28,2.44,2.80,3.44,4.08,4.72,5.36,6.0/ 

c... Minimum vertical size for severe turbulence layer, km 
Data LSEVZ/0. 2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.32, 0.32, 0.32, 

& 0.32,0.25,0.20,0.28,0.36,0.44,0.52,0.60,0.68,0.74,0.82, 

& 0.90,0.98,1.06,1.14,1.22,1.40,1.72,2.04,2.36,2.68,3.0/ 

c... Minimum horizontal size for moderate turbulence layer, km 
Data LMODX/91. ,74. ,60. ,59. ,66. ,69. ,66. ,62. ,54. ,38. ,27. , 

& 22. ,21. ,27. ,35. ,43. ,50. ,58., 66. ,74. ,88. ,92. ,102. ,111. , 

& 121. ,132. ,153. ,200. ,232. ,270. ,300. ,300./ 

c... Minimum horizontal size for severe turbulence layer, km 
Data LSEVX/46. ,37. ,30. ,30. ,33. ,35. ,33. ,31. ,27. ,19. ,14. , 

& 11., 11., 14., 18., 22., 25., 29., 33., 37., 44., 46., 51., 56., 

& 61. ,66. ,77. ,100. ,116. ,135. ,150. ,150./ 


if (model . gt . 7 . or .model . It . 0)stop ’ Invalid Model Number!’ 
c — Non-turbulent case 
if (model . eq. 4) then 
SIGMAU = 0. 

SIGMAW = 0. 

XLWIND = 99999. 

ZLWIND = 99999. 

ISEV = 1 
return 
endif 

sqr2pi = dsqrt ( two/pi ) 
c. . . Find height index for interpolation 

j = 0 

do 10 i = 1,32 

if (Z. lt.hgt(33-i) )goto 10 

j = 33-i 

goto 20 

10 continue 

20 if ( j. It. 1) j = 1 

if( j.gt.31) j = 31 
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c... Interpolate parameters on height 

delz = (z - hgt( j) )/ (hgt( j+1) - hgt(j)) 
c... Average sigmas for moderate turbulence 

smxb = sigmxbar(j) + delz*(sigmxbar( j+l)-sigmxbar( j) ) 
smzb = sigmzbar(j) + delz*(sigmzbar( j+l)-sigmzbar( j) ) 
c. . . Average sigmas for severe turbulence 

ssxb = sigsxbar(j) + delz*(sigsxbar( j+l)-sigsxbar( j) ) 
sszb = sigszbar(j) + delz*(sigszbar( j+l)-sigszbar( j) ) 
c... Scales for sigmas 

XLSIG = lsigx(j) + delz*(lsigx( j+l)-lsigx( j) ) 

ZLSIG = lsigz(j) + delz*(lsigz( j+l)-lsigz( j) ) 
c... Scales for turbulent winds 

XLWIND = lwindx(j) + delz*( lwindx( j+l)-lwindx( j) ) 

ZLWIND = lwindz(j) + delz*(lwindz( j+l)-lwindz( j) ) 
c... Moderate turbulence case 
if (model . eq. 5) then 
SIGMAU = smxb 
SIGMAW = smzb 
ISEV = 2 
return 
endif 

c... Severe turbulence case 
if (model. eq. 6) then 
SIGMAU = ssxb 
SIGMAW = sszb 
ISEV = 3 
return 
endif 

c. . . Interpolate probabilities for encountering turbulence 
pmz = PM(j) + delz*(PM( j+1) - PM(j)) 
psz = PS( j) + delz*(PS( j+1) - PS( j) ) 
c. . . Composite turbulence case 
if (model . eq. 7) then 

SIGMAU = pmz*smxb + psz*ssxb 
SIGMAW = pmz*smzb + psz*sszb 
ISEV = 0 
return 
endif 

c... Minimum thickness for moderate and severe turbulence layers 
maxz(2) = LMODZ(j) + delz*(LMODZ( j+1) - LMODZ ( j ) ) 

maxz(3) = LSEVZ(j) + delz*(LSEVZ( j+1) - LSEVZ(j)) 

maxx ( 2 ) = LMODX(j) + delz*(LMODX( j+1) - LMODX(j)) 

maxx(3) = LSEVX(j) + delz* (LSEVX( j+1) - LSEVX(j)) 


c Rayleigh distribution standard deviations from mean values 

srax = sqr2pi*smxb 
smz = sqr2pi*smzb 
ssx = sqr2pi*ssxb 
ssz = sqr2pi*sszb 

c. . . all = probability of non-turbulent severity category 
all = one - (pmz+psz) 

c. . . als = probability of less than severe category 
als = one - psz 

c. . . Use dx = 0 and dz = 0 to initialize values 
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if (abs(dx) . le. 0. 0. and. abs(dz) . le. 0. 0) then 

c Randomize initial sigma values 

sigmxu = rand(0)*smx 
sigmyu = rand(0)*smx 
sigmxw = rand(0)*smz 
sigmyw - rand(0)*smz 
sigsxu = rand(0)*ssx 
sigsyu = rand(0)*ssx 
sigsxw = rand(0)*ssz 
sigsyw = rand(0)*ssz 
c... Calculate sigma values 
sigu(l) = 0. 
sigw(l) = 0. 

sigu(2) = dsqrt (sigmxu**2 + sigmyu**2) 
sigu(3) = dsqrt (sigsxu**2 + sigsyu**2) 
sigw(2) = dsqrt ( s igmxw**2 + sigmyw**2) 
sigw(3) = dsqrt (sigsxw**2 + sigsyw**2) 
c. . . Moderate or severe turbulence cases 

i f ( mode 1 . eq . 1 . or . mode 1 . eq . 2 ) then 
SIGMAU = sigu(model+l) 

SIGMAW = sigw(model+l) 

ISEV = model+1 
return 
endif 

c. . . Composite turbulence case 
if (model. eq. 3) then 

SIGMAU = pmx*sigu(2) + psx*sigu(3) 

SIGMAW = pmx*sigw(2) + psx*sigw(3) 

ISEV = 0 
return 
endif 

c. . . Select sigma values according to severity level (alpha) 
alpha n rand(0) 

if(alpha. gt. all. and. alpha. It. nls)jsev = 2 
if (alpha. le, all) jsev = 1 
if (alpha. ge. als)jsev = 3 
SIGMAU = sigu(jsev) 

SIGMAW = sigw(jsev) 

c Store turbulence parameter values for next cycle 

psigu = sigmau 
psigw = sigmaw 
ISEV = jsev 
xsev = 0. 
nsev = 1 
return 
endif 

c... Correlations for sigmas 

delx = Dsqrt ( (DX/XLSIG)**2 + (DZ/ZLSIG)**2) 
if (delx. It.0.05d0) then 

rhosig = one - AFAC*delx**2 

else 

rhosig = Dexp(-BFAC*delx) 
endif 

betasig = dsqrt (one - rhosig**2) 
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III. DIGITAL FILTER TECHNIQUES AND TURBULENCE MODELING 


A. BACKGROUND 


The currently used turbulence model for the space shuttle reentry simulation is overly con- 
servative in that severe turbulence is assumed from reentry altitude to the Earth’s surface. In the 
real atmosphere, turbulence is intermittent with large quiescent zones. From reentry to 10 km alti- 
tude, shuttle control is by reaction engines. Johnson Space Center (JSC) sets reaction control 
system redlines based on the predictions of shuttle reentry simulations. Because of the overly con- 
servative model, the orbiter usually lands with about 270 kg (600 lb) of extra fuel. A realistic 
turbulence simulation model was developed to help in more rational selection of redlines. The 
model development was performed in two parts: a stochastic turbulence intensity model and revi- 
sion of the turbulence simulation equations. This paper documents the second part of the develop- 
ment. Part one was performed by another contractor. The two parts are easily integrated into a 
single subroutine; a recommended subroutine structure is presented. 


B. INTRODUCTION 


The current turbulence model for the space shuttle reentry simulation assumes severe turbu- 
lence from reentry to the Earth’s surface. This model is overly conservative because in the real 
atmosphere, turbulence occurs in patches embedded in large quiescent zones. Reaction engines 
provide control of the orbiter from reentry to about 10 km altitude. Fuel budgets for the orbiter are 
based on the reentry simulations with the conservative turbulence model. As a result, the orbiter 
lands with excess fuel. Figure 1 shows the fuel budgets and usage for a typical flight. The orbiter 
takes off with about 2,270 kg (5,000 lb) of fuel. About 45 kg (100 lb) of fuel is used for step- 
away from the external tank after separation. The amount of fuel used in on-orbit maneuvers varies 
and is greatest for rendezvous missions. Currently, on-orbit operations cease when fuel levels drop 
to the 770-kg (1,700-lb) redline. Typically, the shuttle lands with 270 kg (600 lb) of extra fuel. 

For some missions, this redline limits operations. For example, the Hubble Space Telescope (HST) 
mission needs at least 180 kg (400 lb) of additional fuel to put the HST into a higher orbit than is 
normally available from the shuttle. 

A more realistic turbulence model will permit more rational selection of reaction control fuel 
redlines. The task of revising the turbulence model was performed in two parts: development of a 
stochastic turbulent intensity model and revision of the turbulence simulation difference equations. 
The first part was done by Justus [1], The second part is the subject of this study. Part one will be 
discussed briefly, and part two in detail. Most derivations for part two are presented in the 
appendix. All difference equations and associated parameters are presented in a form appropriate 
for coding into flight simulation models. Each equation was tested and verified for accuracy and 
stability. 
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C. TURBULENCE SIMULATION EQUATIONS 


Turbulence is normally simulated as is shown in Figure 2. Gaussian white noise is input to 
the low pass filter. The output simulated turbulence is also Gaussian with the desired spectrum, 
e.g., von Karman. The transfer function of the filter is selected so that the desired output spectrum 
is obtained. If the filter is linear and the probability distribution of the input noise is Gaussian, 
then the output time series will also have a Gaussian probability distribution. Atmospheric turbu- 
lence is often not Gaussian. The linear filter output can be randomly modulated to obtain a more 
realistic output probability distribution. This study deals with generation of the Gaussian turbulence 
rather than with nonlinear modulation of the Gaussian time series. These refinements can be added 
if deemed necessary. 

The most consistently observed characteristic of atmospheric turbulence is the power spectral 
density fall-off with frequency to the -5/3 power. Only occasionally will different fall-offs be 
observed in turbulence spectra. The -5/3 fall-off is consistent with Kolmogorov’s local isotropy 
hypothesis and with von Karman’s spectrum. The difficulty for simulators is the irrational form of 
the von Karman spectrum. In theory, an irrational spectrum gives rise to either an infinite order 
differential equation or a finite, noninteger order differential equation. Investigators have defined 
derivatives of noninteger order and Tatom [2] has used solutions of noninteger order differential 
equations to simulate irrational processes. The computational efficiency of these solutions has not 
yet been demonstrated. 

Frequently, the rational Dryden spectrum, which falls off as frequency to the -2 power, is 
used to simulate turbulence. Frost and Wang [3] have shown quantitative differences between 
simulations of fixed stick aircraft landings flying in Dryden and in von Karman turbulence. The 
differences are in the standard deviation of the touch-down point and in the standard deviation of 
aircraft sink rate. The von Karman spectrum gave lower values of these standard deviations than 
the Dryden model. The Dryden spectrum is, in fact, a crude approximation to the von Karman 
spectrum. The Frost and Wang study provides justification for looking at better approximations. 

Campbell [4] approximated the von Karman spectrum with a higher order rational approxi- 
mation. The discretized Campbell model was shown to have instabilities over some ranges of 
sampling rates. In theory, the Campbell equations were stable over all sampling rates. The 
instability apparently arose from finite precision computer arithmetic. This hypothesis was sup- 
ported by the fact that using extended precision arithmetic delayed the onset of instability. 

Mantey [5] observed a similar phenomenon and pointed out that if the difference equations 
were cast in modal form, the sensitivity of the equations to finite precision arithmetic instabilities 
was minimized. Consequently, a revision of the Campbell formulation was performed that incor- 
porated the Mantey approach. Casting the equations in modal form permitted the analytical deter- 
mination of the dependence of output turbulence standard deviation on time step. Campbell used a 
least squares fit of the observed standard deviation-time step curve in his original paper. 

In this study, a class of transfer function approximations to the von Karman transfer func- 
tion was obtained. For longitudinal spectra, approximations up to fifth order were derived. For the 
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remainder of this paper, the order of an approximation to the von Karman spectrum refers to the 
order of the approximation to the longitudinal spectrum unless the order of the transverse spectrum 
is specifically called out. The corresponding transverse transfer function will be one higher than the 
longitudinal transfer function. 

The dimensionless transfer functions and corresponding spectra are presented in Figure 3. 

The dimensionless transfer functions can be converted to dimensional form by replacing s by aLs/V 
where L is the appropriate turbulent length scale, V is the velocity of the vehicle relative to the 
wind, and “a” is von Karman’s constant (1.339). From the figure, the fifth order approximation 
falls close to the von Karman spectrum. To better understand the closeness of the approximation, 
five simulations were run with the same input Gaussian white noise source data. The fifth order 
approximation was plotted versus the first through the fourth order approximation time series. 

These graphs are shown on the right in Figure 3. If the agreement was perfect, these curves would 
fall along a straight 45-degree line. The greater the scatter about the 45-degree line, the worse the 
agreement. The figure shows that the agreement becomes better as the order of the approximation 
increases. Similar plots were done for different sampling rates. These graphs are not shown, but 
the agreement becomes worse for a given approximation as the sampling rate increases. 

Simulations were performed at several dimensionless sampling rates, with each of the 
approximations and their corresponding transverse equations, to test foi instabilities. These approxi- 
mations were formulated using the Advanced Continuous Simulation Language (ACSL) in single 
precision arithmetic. The results, presented in Tables 1 through 3, show that all simulations are 
stable except for the highest order and the lowest sampling rates. Second and fourth order Runge- 
Kutta and Euler integration were used in the simulations. Since divergence occurred for exception- 
ally long sampling intervals, use of even the fifth order approximation should be feasible for 
shuttle reentry simulations. In any case, the third order simulation was always stable and will be 
presented in detail. 

The difference equations for the third order approximation are presented in Figure 4. The 
equations for both the longitudinal and transverse equations are presented. Longitudinal equations 
will be applied for the along wind, x-component of turbulence. The transverse equations will be 
used for both the y- and z-components. The equations, as shown, can be directly coded into the 
shuttle reentry simulation or into flight simulations for any other vehicle. Required inputs to the 
equations are the appropriate turbulence length scale, the appropriate relative wind velocity, turbu- 
lent intensity, and the sampling rate. The derivation of the equations is presented in the appendix. 
The discretization of the continuous equations was ideal in that output discrete points appear as 
points from the ideally sampled continuous signal. The ideal discretization is easy because of the 
modal form of the equations. The simulation equations currently being used at JSC do not use this 
ideal sampling because of the requirement for speed in a real-time simulation. 
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D. INTEGRATION OF THE TURBULENCE MODEL 


Justus developed a model for turbulent intensities that realistically selects turbulent 
intensities from a Rayleigh distribution. He provided a subroutine that offered eight options for 
turbulent intensities. These options range from zero turbulence to stochastic selection of moderate 
or severe intensities to severe turbulence at all altitudes. For the stochastic case, intensities are 
selected from Rayleigh distributions. Between adjacent altitude bands, the intensities are correlated. 
The Justus model provides all the inputs required by the turbulence equations except for the 
sampling rate. His model has a number of required inputs not currently used in the turbulence 
simulation subroutine used by JSC. The inputs required by the Justus model are shown in Figure 5. 
Justus developed a subroutine that calculates turbulence parameters for the x- and z-components of 
turbulence. The y-component is not explicitly calculated but is generated identically to the z- 
component. The y-component must have its own random number generator but otherwise can use 
coding analogous to that for the z-component. Inputs required for Justus’ subroutine and not cur- 
rently supplied by the JSC model are the changes in position of the vehicle from the last time step 
and the previous values of turbulent intensity. All of these new inputs are readily available from 
the current JSC model. The outputs from the model are the turbulence length scales and standard 
deviations required by the simulation equations. The Justus submodule and the turbulence update 
equations described in the previous section are essentially independent except for the Justus model 
outputs (length scales and standard deviations). As a result, the two efforts are easily integrated. 
The flow chart shown in Figure 6 describes the integration of the two elements into a revised 
subroutine for the JSC simulation. The turbulence equation options available are the original first 
order equations and the third order equations presented here. 


E. RECOMMENDED INVESTIGATIONS 


A higher order, stable approximation to the von Karman spectrum and corresponding differ- 
ence equations were presented in this report. The third order approximation permits simulations 
much closer to real turbulence than simulations based on the rational Dryden spectrum. The 
remaining issue is the relative importance of the higher order approximation to shuttle reentry 
simulation. Simulations based on the third order approximation require much more computation 
than first order simulations. Frost and Wang [3] give some support to the idea that better approxi- 
mations are required. Running a number of simulations with the current shuttle reentry simulation 
using both the current turbulence model and the model with the revised turbulence equations will 
resolve the importance of higher order turbulence models. 

The measure of importance of the higher order simulation will be the relative variability of 
fuel usage for the two turbulence equations. The mean fuel weight after landing will likely be 
about the same for both sets of simulation equations. The standard deviation about the mean will 
be different for the two models. The variation of fuel usage about the mean for the same mission 
affects the on-orbit redline selection. Each set of equations can be used to generate a probability 
distribution of fuel remaining at landing. If the standard deviations of the probability distributions 
are significantly different, use of the results from the higher order simulations are recommended for 
redline selection. 
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x n + 1 
Yn + 1 
WHERE 



DIFFERENCE EQUATION FORM 

-Gx n + Hu n 
■ Cx n 

x„ - STATE VECTOR 
y n > SIMULATED TURBULENT VELOCITY 
u n > GAUSSIAN WHITE NOISE INPUT 
G - COEFFICIENT MATRIX DEFINED BELOW 
H - COLUMN VECTOR DEFINED BELOW 
C - ROW VECTOR DEFINED BELOW 


) 


INPUTS 

o 0 > LONGITUDINAL TURBULENCE STANDARD DEVIATION 
L 0 > LONGITUDINAL TURBULENT LENGTH SCALE 
V > VEHICLE VELOCITY RELATIVE TO THE WIND 
t ■ SAMPLING INTERVAL 


PARAMETERS AND CONSTANTS 


8 ■ 1.339 Lu/V 
X 1 >-1.02025059165 
X 2 >-1.6766157076134 
X 3 >-8.1031337007321 

Ej >exp(XjT/B) 

X. X.2 


c c, 


0.2239020885 

•0.6664166544 

1.442514566 


^3 


1.751871657754 

12.0128342246 

13.86096256684 


Jy Ci2h c 2(Ei-1)2 + 2V V CiCj h Cj h c (Ej-l)(Ej-l)l* 

[6 hnttn M ft XiXj ' (i-e-ej)- . 



MATRIX AND VECTOR DEFINITION 


H 


OuB 

o„ 


h tl (E r 1)/Xi‘ 

hc 2 (E 2 -1)/X 2 

Lh tj (E 3 -1)/X 3 . 




Figure 4a. Equations for the third order longitudinal turbulence model. 
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DIFFERENCE EQUATION FORM 

x„ + i -Gx n + Hu n 

Yn + i -Cx„ 


WHERE 



x n - STATE 

y„ ■ SIMULATED TURBULENT VELOCITY 
u n - GAUSSIAN WHITE NOISE INPUT 
G - COEFFICIENT MATRIX DEFINED BELOW 
H • COLUMN VECTOR DEFINED BELOW 
C - ROW VECTOR DEFINED BELOW 


C INPUTS 

o w - TRANSVERSE TURBULENCE STANDARD DEVIATION 
L w - TRANSVERSE TURBULENT LENGTH SCALE 
V - VEHICLE VELOCITY RELATIVE TO THE WIND 
t - SAMPLING INTERVAL 



PARAMETERS AND CONSTANTS 


B - 1.339 V/L w 
Xi,X 2 , AND X 3 ARE SAME 
AS LONGITUDINAL CASE 
X4 ■•1. 

Ej »exp(XjT/8) 

Cc 2 + _S_ 
X. X.2 


>11.2804723 
>-1.651343027 
> 1 .645595998 
>-10.27472527 


C c + — ♦ — -2- + — 

1 ' ' * X.3 


<<i 

C <2 

C <4 


2.86079443756 

21.368747801 

34.647691313 

13.860962566845 


a _ ^ Cj2 h c 2(Ej-1)2 + 2 V V c l c i h «i h <i (F|*1)(Ej-1) 2 

y X,* (1-6,1) ia fa £ X"xy Tl’-Mj) . 


MATRIX AND VECTOR DEFINITION 

Ei O O O 
O E 2 O O 
O O E 3 O 
O O O E 4 


H 


o wB 

Ou 


hc,(Et-lVXi 

h Cj (E 2 -1)/X 2 

h c ,(E 3 -1)/X 3 

|_h< 4 (E4*1)/X4j 


c - 


1< c 2# c 3» c 4 


Figure 4b. Equations for the third order transverse turbulence model. 




Figure 5. Characteristics of the Justus [1] turbulence model. 
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Figure 6. Recommended turbulence simulation subroutine structure. 





TABLE 1 . STABILITY OF EULER INTEGRATION TURBULENCE SIMULATION 


TRANSFER 

FUNCTION 

ORDER 

ALONG-WIND COMPONENT 
DIMENSIONLESS SAMPLING RATE 

.5 

.1 

.01 

.001 

1 

OK 

OK 

OK 

OK 

2 

OK 

OK 

OK 

OK 

3 

OK 

OK 

OK 

OK 

4 

OK 

OK 

OK 

OK 

5 

DIV* 

OK 

OK 

OK 


TRANVERSE COMPONENT | 

2 

OK 

OK 

OK 

OK 

3 

OK 

OK 

OK 

OK 

4 

OK 

OK 

OK 

OK 

5 

OK 

OK 

OK 

OK 

6 

DIV* 

OK 

OK 

OK 


R4-8255G-104 


* THESE SIMULATIONS DIVERGED 


TABLE 2. STABILITY OF SECOND ORDER RUNGE KUTTA SIMULATIONS 


TRANSFER 

FUNCTION 

ORDER 

ALONG-WIND COMPONENT 
DIMENSIONLESS SAMPLING RATE 

.5 

.1 

.01 

.001 

1 

OK 

OK 

OK 

OK 

2 

OK 

OK 

OK 

OK 

3 

OK 

OK 

OK 

OK 

4 

OK 

OK 

OK 

OK 

$ 

DIV* 

OK 

OK 

OK 


TRANSVERSE COMPONENT f 

2 

OK 

OK 

OK 

OK 

3 

OK 

OK 

OK 

OK 

4 

OK 

OK 

OK 

OK 

S 

OK 

OK 

OK 

OK 

6 

DIV* 

OK 

OK 

OK 


R4-«25SG-105 

• THESE SIMULATIONS DIVERGED 
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APPENDIX 


DERIVATION OF SIMULATION EQUATIONS FOR HIGHER 
ORDER APPROXIMATIONS 


The turbulence simulation equations for the third order approximation in normal and state- 
space form are: 


y 


ttr 


+ k,y" 


+ k 2 y' + k 3 y 


= riu" + r 2 u' + r 3 u 


( 1 ) 


X|' 


-ki 

-k 2 

-k 3 


X| 


1 

X 2 ' 

= 

1 

0 

0 


X 2 

+ 

0 

_X 3 '_ 


0 

1 

0 


X 3 _ 


_ 0 _ 


( 2 ) 


y =? [r t ,r 2 ,r 3 ] 




( 3 ) 


These equations have distinct, real eigenvalues and consequently can be transformed to 
modal form by using the matrix whose columns are eigenvectors of the coefficient matrix. For a 
matrix in the above form, the appropriate transformation is given by: 


T = 


1/X, 

1/X, 2 


1 

l/X 2 

l/\ 2 2 


1 

l/x 3 

1/X 3 2 


( 4 ) 


where X, = eigenvalues. 

The resulting equations are in the form x = Ax + Bu where A is the following diagonal 

matrix. 
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A = 


(5) 


X, 0 0 

0 - A. 2 0 

0 0 \ 3 


The matrix B is given by B = T~'U where U = [1 0 0] T . The superscript T indicates the column 
vector transpose of the row vector [1 0 0], The eigenvalues in equation (5) are in dimensionless 
form. To convert them to dimensional form, the eigenvalues must be divided by 1.339 L/V, L is 
the turbulent length scale and V is the vehicle velocity relative to the wind. To perform a simula- 
tion, the above equations must be discretized. The discretization chosen is in the following form. 


x n+ , = Gx n + Hfu n , 


where 


f = a factor to ensure the correct standard deviation of y 
G = e AT 
H = A -l (e A, -I)B 


t = sampling interval, and 
I = 3x3 identity matrix. 

Once the difference equation is obtained, the variance of the output y is required. For 
equations in the x = Gx + Hu form, the covariance of the state vector is given by the expected 
value of x T x = P. P is given by the solution to the Lyapunov equation 


P = GPG t + HH t 


The unknown quantity in this equation is a matrix. Because G is a diagonal matrix, P can be 
determined easily. 


Pii 


hjhj 

, _ e MP e x .i T/ P 
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where hj is the ith component of H. 

Once P is known, the covariance of y can be determined from: 


CTy 2 = CPC T . 


The preceding equation is expanded in Figure 4. The above approach can be applied for the 
equations of any order. The general equation is: 

a 2 = 2 c i 2 h ci 2 (E i - l) 2 + 2 " n -> CiCiK-jh^Ej - l)(Ej- 1) 

i=| Xi 2 d ~ Ei 2 ) j=7+l j=| XiXjd-EiEj) 

All parameters in the above equation are presented in Figure 4. 

Derivations in this appendix were specifically for the third order equations; the higher order 
equations are in identical form. Only the eigenvalues and the constants are different for different 
equations. 
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